#include <cmath>
#include <vector>
#include <numeric>
#include <algorithm>
#include <functional>
#include "hy1c_util.h"

#ifdef __cplusplus
extern "C" {
#endif
#include "l12_parms.h"
#ifdef __cplusplus
}
#endif

namespace hy1c {
    template <typename T>
    struct Type
    {
        typedef T (*fun)(T a[], size_t n);
    };


    template <typename T>
    void swap(T &a, T &b)
    {
        T t = a;
        a = b;
        b = t;
    }

    template <typename T>
    void partition(T A[], size_t left, size_t right, size_t &pos)
    {
        if (A == nullptr)
            return;

        T data = A[left];

        size_t i = 0;
        for (pos = left, i = left + 1; i <= right; i++)
        {
            if (A[i] < data)
            {
                pos++;
                swap(A[pos], A[i]);
            }
        }

        swap(A[left], A[pos]);
    }

    template <typename T>
    T getmid(T A[], size_t n)
    {
        if (n <= 0)
            return 0;

        size_t left = 0;
        size_t right = n - 1;
        size_t mid = (left + right) / 2;
        size_t pos = 0;

        while(true)
        {
            partition(A, left, right, pos);
            if (pos == mid)
                break;
            else if (pos > mid)
                right = pos - 1;
            else
                left = pos + 1;
            //1th  pos=3   1 5 3 7 11 9
            //2th  pos=0   1 5 3 7 11 9
            //3th  pos=2   1 3 5 7 11 9
        }

        return (n & 0x1) != 0 ? A[mid] : (A[mid] + A[mid + 1]) / 2;
    }

    template <typename T>
    T average(T A[], size_t n)
    {
        T sum = std::accumulate(A, A + n, 0.0);
        return sum / n;
    }

    // standard deviation
    template <typename T, typename F>
    T stdev(T a[], size_t n, F fun = average)
    {
        T mean = fun(a, n);

        T accum;
        std::for_each(a, a + n, [&](const T &d) {
                accum += pow((d - mean), 2);
                });

        T stdev = sqrt(accum / n);
        return stdev;
    }

    template <typename T>
    int line_interpolation(T * wave, size_t nwave,
            T * wave_ref, size_t nwave_ref,
            T * in, T * out)
    {
        for (size_t i = 0; i < nwave; i++) {
            T x1 = wave_ref[0];
            T x2 = wave_ref[1];
            T y1 = in[0];
            T y2 = in[1];

            int wave_equal = 0;

            for (size_t j = 0; j < nwave_ref; j++) {
                if (wave[i] == wave_ref[j]) {
                    wave_equal = 1;
                    out[i] = in[j];
                    break;
                } else if ((wave[i] < wave_ref[j]) || (j == nwave_ref - 1)) {
                    if (j > 1) {
                        x1 = wave_ref[j - 1];
                        x2 = wave_ref[j];
                        y1 = in[j - 1];
                        y2 = in[j];
                    }
                    break;
                }
            }

            if (wave_equal)
                continue;

            T a = (y2 - y1) / (x2 - x1);
            T b = y2 - a * x2;
            out[i] = a * wave[i] + b;

            //printf("in %f out %f\n", in[i], out[i]);
        }
        return 0;
    }

    template <typename T>
    T frequency_distribution_histogram(T * a, size_t n)
    {
        T max = *std::max_element(a, a + n) * 1.01;
        T min = *std::min_element(a, a + n) * 0.99;
        constexpr size_t slice_cnt = 5; // TODO configuration
        const float offset = (max - min) / slice_cnt;

        //size_t his_array[slice_cnt] = {0};
        std::vector<std::vector<T> > his_array(slice_cnt);

        for (int i = 0; i < n; i++) {
            int index = static_cast<size_t>((a[i] - min) / offset);
            his_array[index].push_back(a[i]);
        }

        for (int i = 0; i < slice_cnt; i++) {
            for (auto j = 0; j < his_array[i].size(); j++)
                printf("#");
            printf(" %d\n", his_array[i].size());
        }

        auto ipeak = *std::max_element(his_array.begin(), his_array.end(),
                [](auto a, auto b) -> bool {
                return a.size() < b.size();
                });

        T peak = std::accumulate(ipeak.begin(), ipeak.end(), 0.0) / ipeak.size();
        //std::for_each(ipeak.begin(), ipeak.end(), [](const T &a){printf("%f ", a);});
        //T peak = offset * (2 * (ipeak - his_array.begin()) + 1) / 2 + min;
        printf("\n%f\n", peak);
        printf("\n");
        return peak;
    }




} // namespace hy1c

float getmid_float(float a[], size_t n)
{
    return hy1c::getmid<float>(a, n);
}

float stdev_float(float a[], size_t n)
{
    return hy1c::stdev<float, hy1c::Type<float>::fun>(a, n, hy1c::getmid);
}

float line_interpolation_float(float * x, size_t n,
        float * x_ref, size_t n_ref,
        float *in, float *out)
{
    return hy1c::line_interpolation<float>(x, n, x_ref, n_ref, in, out);
}

float frequency_distribution_histogram_float(float * a, size_t n)
{
    return hy1c::frequency_distribution_histogram<float>(a, n);
}

#if 0
#include <iostream>

int main(int argc, char *argv[], char *env[])
{
    float Array[] = {7.5, 5.3, 3, 11, 9, 3.14, 3.14, 3.14, 3.14};

    hy1c::frequency_distribution_histogram<float>(Array, sizeof(Array)/sizeof(Array[0]));
    //	float mid;
    //	mid = getmid_float(Array, sizeof(Array) / sizeof(Array[0]));
    //
    //	float stdev = hy1c::stdev<float>(Array, sizeof(Array) / sizeof(Array[0]));
    //	std::cout << "The mid number of the array is: " << mid << std::endl;
    //	std::cout << "The stdev of the array is: " << stdev << std::endl;

    return 0;
}
#endif
